How dense can one pack spheres of arbitrary size distribution? 
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PACS 81.05.Rm - Porous materials; granular materials 

PACS 45 . 70 . -n - Granular systems 

PACS 45.70.Cc - Static sandpiles; granular compaction 

Abstract - We present the first systematic algorithm to estimate the maximum packing density 
of spheres when the grain sizes are drawn from an arbitrary size distribution. With an Apollonian 
hlling rule, we implement our technique for disks in 2d and spheres in 3d. As expected, the densest 
packing is achieved with power-law size distributions. We also test the method on homogeneous 
and on empirical real distributions, and we propose a scheme to obtain experimentally accessible 
distributions of grain sizes with low porosity. Our method should be helpful in the development 
of ultra-strong ceramics and high performance concrete. 



High strength ceramics and ultra-high performance con- 
crete (UHPC) [l] require minimizing the porosity out of 
a compacted either reactive or sintered powder. Much re- 
search effort has been invested in the past to optimize this 
procedure and the most important factor turned out to 
lie in the adequate choice of the size distribution of the 
constituents. In fact, the broader this distribution, the 
smaller result the remaining voids. So, mixtures of very 
different grains with up to four orders of magnitude in 
grain size are typically used for UHPC. But what is the 
ideal partitioning? Which combined grain-size distribu- 
tion would yield the most compact mixture? The key to 
answer this fundamental question posed by the practition- 
ers is to be able to estimate the maximum density a given 
distribution can provide. This is precisely the aim of this 
Letter. 

The maximal filling density has been studied for many 
specific types of packings [2] starting with the work of 
Fuller [3] . Besides the simple monodisperse case also exact 
results for some symmetric cases are known ^|5| . Various 
models have been proposed to deal with the superposi- 
tion of two or three rather sharply peaked distributions, 
like the ones by Toufar et al. [6], Yu and Standish |7j, 
or the various linear theories by De Larrard [8|[9]. Also 
for continuous size distributions a hierarchical partition- 
ing model was recently developed [To]. Most real grain- 
size distributions, used for dense packing, have a rather 
complicated shape often being a superposition of various 



empirical functions. It is therefore of interest to develop a 
technique to obtain an upper bound for a packing having 
arbitrary distribution. 

Based on the insight that a completely space-filling 
packing of spheres, i.e., having unity volume density, can 
only be achieved with a generalized (random) Apollonian 
setup 11 12 , we design a systematic technique to op- 
timally fill the fines into the voids between larger grains 
and, by sweeping from the large end of any distribution, to 
finally obtain the highest density one could possibly attain 
with such distribution. Testing the technique on various 
real and artificial distributions we recover that power laws 
provide the highest densities. 

As illustrated in Fig. [T] we discretize the distribution 
into bins by grain size. Grains are assumed to have spher- 
ical shape (disks in two dimensions) and are organized 
in each bin in the monodisperse closest-packing configu- 
ration. Then the gaps are filled with smaller ones, fol- 
lowing an Apollonian packing. As rigorously proved by 
L. F. Toth [13 14 , in two dimensions, the most efficient 
monodisperse arrangement of disks is the hexagonal closest 
packing [hep), with a density of phcp = l/67r-\/3 « 0.9069. 
With this arrangement, there are, per largest grain, two 
unitary voids to be filled with the traditional Apollonian 
packing (see Fig. [T]). Each bin b is characterized by four 
different parameters: the average radius, r;,; the volume 
of grains, VJ,; the effective volume, V^f^\ and the density 
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of the arrangement of particles in the bin, pi,. Since we 
consider all grains to have the same density the volume 
of grains, Vf,, quantifies the total amount of matter. The 
excluded volume interaction between grains imposes ge- 
ometrical restrictions such that, even for an efficient ar- 
rangement, the effective volume, VJ,g, is larger than VJ,, 
accounting for the volume of both the grains and voids. 
The density of the arrangement in the bin, pf,, quantifies 
the efficiency of the packing and can be related with the 
volumes as Vb = PbV^g. 

Figure [T] is a pictorial scheme of the sequential com- 
pacting procedure, with bins representing grains grouped 
according to their size. The height of each bin represents 
the effective volume, V^*^, and bins are organized in the 
inverse order of their radius. Initially, all bins are consid- 
ered to be in the closest-packing configuration, hep, i.e., 
their density is pb = Phcp- Starting with the bin 1, corre- 
sponding to the largest grains, both the radius and volume 
of grains from further bins, required to fill the voids, are 
computed and transferred to the bin (details below). The 
volume of grains, Vi, and the density, pi, are updated to- 
gether with the effective volume and volume of grains of 
the bins from which particles are transferred. The process 
is executed for each bin, from the largest to the smallest 
grains, and the final net density, Pnot, is then computed as 



Pnet 



(1) 



where Vf, is the sum of the initial volume of particles in 
the bin with the volume of particles collected from the fol- 
lowing ones. In the denominator, V^g is the final effective 
volume of the bin. 



In the two dimensional case, for each generation of the 
Apollonian packing, the radius and the number of grains 
(disks) collected from other bins can be obtained, from the 



positive solutions of the Soddy-Gossett equation 15 



(2) 



where the sum runs over the three disks limiting the void 
having radius, r,;, and rj is the radius of the required one. 
The necessary volume of disks, needed to fill all voids, is 
then transferred from the corresponding bin. If not enough 
volume is available, we proceed to the following bins until 
the required number of grains is collected. 

Let us consider, as example, a grain-size distribution 
such that the volume of grains, Vb, is a power law of the 
grain size, i.e., Vb ^ (l/r)~". As in real situations, the 
grain-size distribution is truncated at a lower cutoff, e. 
Once initially grouped in bins, by size, and arranged in 
hep, Vcs ^ (l/r)~". Figure [2] shows the initial (circles) 
and final (full line) effective volume, VeS, as a function of 



the inverse radius, for a — 0.71 and e ~ 10^^. The pro- 
posed algorithm gives a final net density pnct ~ 0.9997, 
which is close to unity and a porosity 300 times smaller 
then the closest-packing density of monodisperse disks. 
The effective volume of smaller-grain bins is significantly 
reduced and, for the tiny grains, the effective volume van- 
ishes, i.e., all material is completely used to fill previous 
bins. The inset shows the final density of each bin. The 
density of the first bins is close to unity and the pack- 
ing efficiency is solely limited by the cutoff. The smaller 
the grain size the lower the density since the number 
of Apollonian-packing generations which can be collected 
from further bins diminishes and the reserve from such 
bins also vanishes. The density of the final bins can either 
correspond to the one of the closest packing or be zero in 
the case that no material is left. In the same inset (Fig. [2]) 
we also show the data for different cutoffs. For all con- 
sidered values the same qualitative behavior is observed 
but the smaller the cutoff (minimum radius) the larger 
the number of different bins with density close to unity 
and the more effective the packing. 

Figure [3] shows the dependence of the final net density 
on the exponent a, for different values of the cutoff, e. 
For each cutoff, we observe an optimal value of a = aopt 
at which a maximum final net density is obtained, that 
decreases with the cutoff, as shown in the inset of Fig. [3] 
In the limit of vanishing cutoff, the optimal exponent con- 
verges to 2 — df, where d/ is the fractal dimension of the 
Apollonian packing, which in 2d is df « 1.306 16 . In 
this optimal case, the density is unity and the distribu- 
tion of grain volume is the one of the Apollonian pack- 
There exist different families of determin- 
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istic Apollonian packing with fractal dimensions ranging 
from 1.306 to 1.802 [ll], which will lead to different val- 
ues Qfopt = d — df {d is the dimension of the system). 
Random Apollonian packings, instead, are characterized 
by the same fractal dimension (as shown by Baram and 
Herrmann '12|) and have, therefore, the same Ofopt- 

In 3d, grains are spheres and, for the sake of conve- 
nience, we start with a regular tetrahedron configuration 
with four mutually tangent spheres that are also tangent 
to an enveloping sphere, with a density p w 0.3633 18 . 
Recently, Baram and Herrmann 19 20 have introduced 



an algorithm to construct an Apollonian packing in 3d 
starting from any initial configuration, which we consider 
here to compute the radius and the number of grains to 
collect from further bins. Alike the 2d case, when an ini- 
tial power-law distribution of volumes is considered, with 
a = 0.55 and cutoff £ = 3 x lO^'^, the final net den- 
sity is p ~ 0.9325, corresponding to a porosity that is 
1/3 of the one for the closest packing of monodisperse 
spheres. The aopt converges to 3 — df, with df « 2.4739 
the fractal dimension of the Apollonian packing 18 . For 
a = dopt , the porosity p scales with the cutoff as p ^ e"^ , 
with T — 0.47 ± 0.03, close to the one of the Apollonian 
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In real appHcations, since a fine control of the grain size 
is not feasible, the usual distributions differ significantly 
from the idealized power laws. Due to the discrete na- 
ture of the system, the optimal distribution would be a 
sequence of delta functions centered at the characteristic 
grain sizes of the Apollonian arrangement and with pref- 
actors given by a power law with exponent aopt. But, such 
extremely narrow distributions are experimentally not re- 
alized and, instead, Gaussians with a certain size disper- 
sion typically appear as the generic distribution in prati- 
cal situations and are used here to properly describe real- 
istic grain-size distributions. Partitioners typically con- 
sider mixtures of materials with different characteristic 
grain sizes like, e.g., crushed quartz, cement, sand, sil- 
ica fume, and microsilica, in the case of high performance 
concrete (HPC). To illustrate such procedure, in Fig.|4]we 
take a representative empirical distribution for HPC, with 
a mixture of four components, obtained from Ref. 21 1, 
for which we estimate an upper bound for the density of 
Pjjpj « 0.8203 which is consistent with the typical val- 
ues discussed in the reference (around 0.8). Optimizing 
the grain-size distribution, by finding an efficient set of 
different sizes to minimize the porosity, is a relevant tech- 
nological problem that can be addressed in a systematic 
way with our algorithm. For each component i, the dis- 
tribution is characterized by a Gaussian with average size 
yUi, size dispersion tJi, and height of the peak Hi. We start 
by considering two types of material, i.e., two Gaussians. 
For simplicity, to reduce the number of degrees of freedom 
we fix their volume fraction to be equal. By exploring 
the parameter space, the algorithm can identify the ra- 
tios between average sizes and size dispersions which min- 
imize the porosity (see Fig. [5^). Once this distribution has 
been identified, a third Gaussian can be included and its 
space of three parameters explored (see Fig. [5|d) . The pro- 
cess can be repeated until the desired porosity is obtained. 
This gives an easy recipe to a systematic proportioning of 
mixtures for highly compacted powders. 



ticeable improvement, whereas in 3d (with a ~ 0.55 and 
£ = 3 X 10^"^), since the packing efficiency is significantly 
affected by these effects, popt ~ 0.9493. The latter corre- 
sponds to a decrease in almost 2% in the porosity. The 
effect is even more significant with the empirical distribu- 
tion in Fig. |4] In this case, a more efficient distribution is 
obtained (represented by the gray area in the plot), with 
a maximum density popt ~ 0.9501, which is three times 
lower in porosity than the original case. The proposed op- 
timization technique could be implemented in practice for 
instance by using adequate filtering. 

In summary, we propose an iterative process to esti- 
mate the upper bound for the density of a packing of 
spheres having an arbitrary size distribution, which is 
meaningful for developing low-porosity materials. Grains 
are grouped by their sizes and, sweeping from the largest 
to the smallest grain size, the voids are hierarchically filled 
with smaller grains according to an Apollonian packing. 
We have analyzed the dependence of the optimal density 
on the properties of the distribution and recovered the 
power laws giving the most efficient case. We have sug- 
gested an iterative scheme to optimize the proportioning of 
the components in order to reach the lowest porosity. Fu- 
ture work should consider different shapes of grains or even 
a broad distribution of shapes. Developing experimental 
hierarchical procedures to obtain the efficient packings re- 
ported would be of paramount interest. 
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Fig. 1: Schematic description of the algorithm in 2d. a) Using 
the hep as the initial packing of a particular bin, the algorithm 
fills the empty spaces using the positive solution of Eq. ([2|. b) 
By symmetry, for each largest disk, two voids are to be filled, 
c) In an arbitrary distribution of effective volumes (right-hand 
side), we remove the amount of volume required according to 
the procedure described in the text. In the beginning all bins 
are filled with the hep. Starting from the bin of the largest 
grains, voids are filled with grains from the further bins: forth 
bin (red) for the first generation, seventh bin (blue) for the 
second, and tenth and eleventh bins (magenta) for the third. 
The regions shown in the first bin represent the final contri- 
bution of each grain size to the effective volume, where white 
stands for the empty space. Proportion between sizes has been 
exaggerated for clarity. 
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Fig. 2: Initial and final distribution of tlie effective volume 
as a function of 1/r, in 2d. The circles correspond to the 
initial distribution, given by Vcs ~ (l/r)~°, with a = 0.71 and 
cutoff £ — 10~^ (minimum radius). The full-(blue) line, is the 
final distribution, with a density of Pnet- The inset shows the 
density of each bin as a function of 1/r (diamonds). Data for 
e = 10~^ (squares) and e = 10~* (triangles) is also included. 
The dashed- (blue) line stands for the minimum density phcp ~ 
0.9069, which corresponds to the initial configuration. 




'I ' I 1 



= 0.3633 

10"' W 10' io" 
7/r[l/iim] 



q <?q 



. (DO 



IlJ I ' I ' 



10 



10^ 10' 

1/r [l/|im] 



10 



Fig. 4: Initial and final distribution of the effective volume as 
a function of 1/r, for 3d. The initial values were obtained from 
the empirical distribution of ultra-high performance concrete 
of Fig. 7 of Ref. |21j . The circles correspond to the initial distri- 
bution and the full- (blue) line is the final distribution, giving a 
density of pnct ~ 0.8203. The gray area stands for the optimal 
distribution proposed in the text, giving popt ~ 0.9501. As in 
Fig. [2] the inset shows the density of each bin as a function 
of 1/r (diamonds). The dashed-(blue) line stands for the min- 
imum density, which corresponds to the initial configuration. 



I^llllllj I I lllll^ I I IMIM| I I llllll 




Fig. 3: Final net density, pnot, as a function of q, for different 
values of the cutoff e, in 2d. All curves are characterized by an 
optimal a = Qopt, at which a maximum density is obtained. 
The inset shows the dependence of ctopt on the inverse of the 
cutoff e. When the effect of truncation is reduced the value of 
«opt approaches 2 — df (dashed line) as Ofopt — (2 — d/) ~ l/e"", 
where d/ is the fractal dimension of the Apollonian packing and 
a = 0.37 ±0.09. 
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Fig. 5: a) Final net density, Pnot, for two Gaussians as a 
function of the ratio between size dispersions ((T2/cri). Each 
set corresponds to different relative average sizes {1x2/ fii)- b) 
Final net density for three Gaussians, where the first two have 
Ma/pi = 3 and (J2/cri = 2.1 (the highest density in a)), as a 
function of the size dispersion (crs/cri) and for different average 
sizes (pa/pi) of the third one. For both cases, e = 3 x 10~^. 
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